浅谈C++浮点数精度的一些问题

有时，浮点数的精度是在做题中的一个重要的问题。举一个例子：

#include <stdio.h>
main() {
    printf("%.6lf\n", floor(sqrt(9999.9999) * sqrt(10000.0001)));
}

这份代码是一个简单的根号的乘法，我们可以手算得到 $\sqrt{9999.9999 \times 10000.0001} = \sqrt{10000^2 - 0.0001^2} < 10000$，应当输出 9999.000000 。但是当你真正运行的时候，会发现它会输出 10000.000000 ，而不是 9999.000000 。这是为什么呢？这就要从浮点数的精度说起。

1. 浮点数精度和范围

在说浮点数精度之前，我们先说说一个相关的东西：科学计数法

用科学计数法表示的数一般是这样子的：$\text{符号} a \times b^c$ ，其中符号为正或负数，$a$ 是有效数字，$b$ 是底数，$c$ 是指数。而且可以很容易知道，$1 \le a < 10$ 。例如在 $-3.14 \times 10^{15}$ 中，- 是符号，3.14 是有效数字，10 是底数，15 是指数。

在计算机中，浮点数也用科学计数法表示。不过，因为用的是二进制存储，所以计算机中的浮点数一般这样表示：

$\text{符号} a \times 2^c$

例如，二进制下的 $10101.101$ （对应十进制下的 $21.625$）就可以用科学计数法表示为$(1.0101101)_2 \times 2^4$

在C++中， float 是单精度，32位，double 是双精度，64位。并且这两个在内存中的分配结构也不一样：

这里以 float 为例，符号位占 $1$ 位，$0$ 表示正数，$1$ 表示负数。

指数占 $8$ 位，表示 $0 \sim 255$ ，为了能表示负指数，实际指数位存储在指数域中值减去一个偏移量（float 为 127 ，double 为 1023），例如用 0000 0000 表示 0-127=127 ，1000 0000 表示 128-127=1 ，1001 0001 表示 145-127=18 ，0001 0001 表示 17-127=-110 ，以此类推。（采用了移位存储）

尾数，也就是科学计数法中的底数，占 $23$ 位。显然底数的整数部分必定是 $1$ ，也就是这个小数必定是 $1. \cdots \times 2^c$ ，所以记的时候就把前面的 $1$ 省略掉了，只记后面的小数部分。

可以发现，由于 $2^{23} = 8388608$ ，一共 $7$ 位，这意味着最多能有 $8$ 位有效数字（还要算上前面的 $1$），但绝对能保证的为 $7$ 位，也就是说 float 的精度为 $7 \sim 8$ 位有效数字；并且表示范围为 $[-2^{128}, 2^{128})$ ，也就是 $-3.40 \times 10^{38} \sim 3.40 \times 10^{38}$

同理，由于 $2^{52} = 4503599627370496$ ，一共 $16$ 位，这意味着最多能有 $17$ 位有效数字（还要算上前面的 $1$），但绝对能保证的为 $16$ 位，也就是说 float 的精度为 $16 \sim 17$ 位有效数字；并且表示范围为 $[-2^{1024}, 2^{1024})$ ，也就是 $-1.80 \times 10^{308} \sim 1.80 \times 10^{308}$

举个例子：$(12.625)_{10} = (1100.101)_2 = (1.100101)_2 \times (2^3)_{10}$ ，所以这个浮点数在双精度下的存储是：

2. 进制转换

整数的进制转换相信大家都会，那么这里只讲小数的进制转换。这里，小数的进制转换你会发现有三种情况：

(1) 有限小数变成有限小数 ，举个最简单的例子：$(12.5)_{10} = (30.2)_{4}$

(2) 有限小数变成无限小数（或反过来） ，举个例子：$(1.2)_{10} = (1.00110011 \cdots)_2$ ，$(\frac{10}{3})_{10} = (10.1)_3$

(3) 无限小数变成无限小数，举个例子：$(\frac{10}{3})_{10} = (10.010101\cdots)_{2}$

通过这些例子，可以发现，有些小数是无限的。

3. 精度失真

有些小数是无限的，但是精度是有限的，就会导致经过处理后的数据并不准确。这里我们举个例子方便理解：

$(2.2)_{10} = (10.00110011 \cdots)_2 = (1.000110011 \cdots)_2 \times (2^1)_{10}$ ，由于精度有限，转换后如下表所示：

而上表中的数转换成十进制之后已经不是 $2.2$ 了，而是 $2.200000000000000200$ 。可以通过以下代码进行测试：

#include <stdio.h>
main() {
    printf("%.18lf\n", 2.2);
}

4. 例子的解释

我们来看一开始举的例子：

#include <stdio.h>
main() {
    printf("%.6lf\n", floor(sqrt(9999.9999) * sqrt(10000.0001)));
}

通过上面的解释，可以发现：

$9999.9999 \rightarrow 9999.9999000000007$

$10000.0001 \rightarrow 10000.000099999999$

开完算数平方根后（左边是失真的数开根后实际的数，右边是失真后的数）：

$99.999999500000002285000011425 \rightarrow 99.9999995000000012$

$100.000000499999995200000024 \rightarrow 100.0000004999999987$

相乘后得到（左边是失真的数相乘后实际的数，右边是失真后的数）：

$9999.9999999999998500 \rightarrow 10000.0$

于是，原本到不了 $10000$ 的数变到了 $10000$ 。

5. 解决方法

其中一种解决的方法就是少用除法和少开根号 。假如你要比较 $\frac{a}{b}$ 和 $\frac{c}{d}$ 的大小，你就可以通过比较 $a \times d$ 和 $c \times b$ 的大小来减小失真（可以使得原来是有限小数的数（二进制下）还是有限小数，失真的概率小些）。假如你要比较 $\sqrt{x^2 + y^2}$ 与 $z \text{ } (z > 0)$ 的大小，你就可以直接比较 $x^2 + y^2$ 和 $z^2$ 的大小来减小失真。

还有一种解决方法就是手写浮点数，像这样：

...
struct BigDecimal {
    int before, after; //  分别表示小数点前的十进制数 和 小数点后的十进制数
    int pm = 0; //表示正负号, 这里 0 为正, 1 为负
    int length = 30, base = 1e8 + 0.5; // 表示小数点的长度 和 小数点后的进位制 (这里是 after 满 10^8 向前进位)
    BigDecimal(int pp = 0, int bb = 0, int aa = 0) {
        pm = pp; before = bb; after = aa;
    }
    bool operator < (const BigDecimal &A) const {
        if (pm != A.pm) return pm > A.pm;
        if (pm == 0) return before != A.before ? before < A.before : after < A.after;
        return before != A.before ? before > A.before : after > A.after;
    }
    bool operator > (const BigDecimal &A) const {
        return A < (*this);
    }
    bool operator == (const BigDecimal &A) const {
        return !((*this) < A || A < (*this));
    }
    ...
};
BigDecimal operator + (const BigDecimal A, const BigDecimal B) {
    if (A.pm == B.pm) return BigDecimal(A.pm, A.before + B.before + (A.after + B.after) / A.base, (A.after + B.after) % A.base);
    else ...
}
...

这里的原理就是本来小数点后的数会失真，将其变为整数后就不会失真了。

或者有的时候这样写可以更加简单（这里假设保留 $4$ 位小数）：

...
typedef long long LL;
const LL D = 1e4 + 0.5;
LL input() {
    double x;
    scanf("%lf", &x);
    return llround(x * D); // llround()  指的是返回值为 long long, 而参数没有变化
}

功能：将这个数的小数点往右移动 $4$ 位（最后一位要四舍五入），将一个小数变为一个整数，减少失真。（假如想要算出这个数，就将这个数的小数点往左移动 $4$ 位）

那你可能就会问了：假如我就是要算出来根号怎么办？这不难，二分啊！（这里以已知直角边算出斜边为例）

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
LL A, B, C;
int main() {
    scanf("%lld%lld", &A, &B);
    LL l = 0, r = A * A + B * B;
    while (l <= r) {
        LL mid = (l + r) >> 1;
        if (mid * mid <= A * A + B * B) l = mid + 1;
        else r = mid - 1;
    }
    printf("%lld\n", r); // r 是向下取整, l 基本上是向上取整 (l = r + 1) , 但在 r * r = A * A + B * B 的时候，l 就不是向上取整了.
    return 0;
}

这里以一道题为例 ABC191D ：给你一个圆的圆心坐标 $(X, Y)$ ，以及圆的半径 $R$ ，问你园内整点的个数。（ $\vert X \vert ,\vert Y \vert \le 10^5$ ，并且 $X, Y$ 小数点后最多有 $4$ 位）

如果用浮点数就会被卡精度，那么我们就这么写：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
LL read() {
    double x;
    scanf("%lf", &x);
    return llround(x * 10000);
}
LL sq(LL x) { return x * x; }
int main() {
    X = read(); Y = read(); R = read();
    LL ans = 0, RR = sq(R);
    LL high = Y+1, low = Y-1, pos = X - R;
    while (pos % 10000) pos++;
    while (high % 10000) high++;
    while (low % 10000) low--;
    for (; pos <= X + R; pos += 10000) {
        LL A = sq(pos - X);
        while (high >= Y && A + sq(high - Y) > RR) high -= 10000;
        while (high <= Y || A + sq(high - Y) <= RR) high += 10000;
        while (low <= Y && A + sq(low - Y) > RR) low += 10000;
        while (low >= Y || A + sq(low - Y) <= RR) low -= 10000;
        LL add = (high - low) / 10000 - 1;
        if (add > 0) {
            ans += add;
        }
    }
    printf("%lld\n", ans);
    return 0;
}